import math
import numpy
#from libpy.model_mddataset import MdDataset
from libpy.model_mdobject import MdObject
from statlib import stats 
from numpy import newaxis as nA
import libpy.chemometrics 

class MdDatamatrix:
    def __init__(self):
        pass
    def AddDataset(self,dataset):
        self.dataset = dataset
        self.dimension = dataset.dimension
        self.nObservation = len( dataset.objects )
        self.nVariable = len( dataset.objects[0].landmarks ) * self.dimension
        self.matrix = numpy.zeros( ( self.nVariable, self.nObservation ) )
        i = 0
        #for object in dataset.objects:
        for obj in dataset.objects:
            j = 0
            #for lm in object.landmarks:
            
            for lm in obj.landmarks:
                #print lm.xcoord, lm.ycoord, lm.zcoord
                self.matrix[j,i] = lm.xcoord
                j += 1
                self.matrix[j,i] = lm.ycoord
                j += 1
                if self.dimension == 3:
                    self.matrix[j,i] = lm.zcoord
                    j += 1
                if j == self.nVariable:
                    break
            i += 1
    
class MdManova:
  def __init__(self):
    self.dimension = -1
    self.data = MdDatamatrix()

  def AddDataset(self,dataset):
    self.data.AddDataset( dataset )
    self.groupinfo_idx = 0
    self.group_list = {}
    i = 0
    for object in self.data.dataset.objects:
      group_key = object.group_list[self.groupinfo_idx]
      print group_key
      if( not self.group_list.has_key(group_key) ):
        self.group_list[group_key] = []
      self.group_list[group_key].append( i )
      i+=1
    return

  def Analyze(self):
    variable_index = 0
    data = {}
    data['sum'] = [0] * self.data.nVariable
    data['count'] = 0
    self.group_data = {}
    for key in self.group_list.keys():
      self.group_data[key] = {}
      within_group_sum = [0] * self.data.nVariable
      for observation_index in self.group_list[key]:
        for variable_index in range( self.data.nVariable ):
          val = self.data.matrix[variable_index,observation_index]
          within_group_sum[variable_index] += val
          data['sum'][variable_index] += val
      self.group_data[key]['sum'] = within_group_sum
      self.group_data[key]['count'] = len( self.group_list[key] )
      self.group_data[key]['average'] = []
      for sum in within_group_sum:
        self.group_data[key]['average'].append( sum / float( self.group_data[key]['count'] ) )
      data['count'] += self.group_data[key]['count']
      #print self.group_data[key]

    data['average'] = []
    for sum in data['sum']:
      data['average'].append( sum / float( data['count'] ) )
    W = numpy.zeros( ( self.data.nVariable, self.data.nVariable) )
    #nGroups = len( self.group_list.keys() )
    B = numpy.zeros( ( self.data.nVariable, self.data.nVariable ) )
    x = []
    group = []
    for key in self.group_list.keys():
      for variable_index_i in range( self.data.nVariable ):
        for variable_index_j in range( self.data.nVariable ):
          val_i = self.group_data[key]['average'][variable_index_i]
          val_j = self.group_data[key]['average'][variable_index_j]
          B[variable_index_i,variable_index_j] += self.group_data[key]['count'] * ( val_i - data['average'][variable_index_i] ) * ( val_j - data['average'][variable_index_j] )
      print B
      for observation_index in self.group_list[key]:
        x.append( [] )
        group.append( key ) 
        for variable_index_i in range( self.data.nVariable ):
          x[-1].append( self.data.matrix[variable_index_i,observation_index])
          for variable_index_j in range( self.data.nVariable ):
            val_i = self.data.matrix[variable_index_i,observation_index]
            val_j = self.data.matrix[variable_index_j,observation_index]
            W[variable_index_i, variable_index_j] += ( val_i - self.group_data[key]['average'][variable_index_i] ) * ( val_j - self.group_data[key]['average'][variable_index_j] )
      print W
    print x
    print group
    b, w = _BW( numpy.array(x), numpy.array(group) )
    print b
    print w
    scores, loads, eigs = libpy.chemometrics.cva( numpy.array(x), numpy.array(group), None)
    print "scores", scores
    print "loads", loads
    print "eigs", eigs
    '''
    # so.. I have between group sums of squares and within group sum of squares. now what? 
    df1 = ( len( self.group_list.keys() ) - 1 )
    bgms = bgss / df1
    print "bgms", bgms
    df2= ( data['count'] - len( self.group_list.keys() ) )
    wgms = wgss / df2
    print "wgms", wgms
    f = bgms/wgms
    print "F:", f
    prob = stats.fprob(df1,df2,f)
    print "p val", prob
    f, prob = stats.lF_oneway( anova_data[0], anova_data[1] )
    print f, prob
    return
    '''

def _index(y,num):
    """use this to get tuple index for take"""
    idx = []
    for i in range(len(y)):
        if y[i] == num:
            idx.append(int(i))
    return tuple(idx)

def _BW(X,group):
    """Generate B and W matrices for CVA
    Ref. Krzanowski
    """
    mx = numpy.mean(X,0)[nA,:]
    tgrp = numpy.unique(group)
    for x in range(len(tgrp)):
        idx = _index(group,tgrp[x])
        L = len(idx)
        meani = numpy.mean(numpy.take(X,idx,0),0)
        meani = numpy.resize(meani,(len(idx),X.shape[1]))
        A = numpy.mean(numpy.take(X,idx,0),0) - mx
        C = numpy.take(X,idx,0) - meani
        if x == 1:
            Bo = L*numpy.dot(numpy.transpose(A),A)
            Wo = numpy.dot(numpy.transpose(C),C)
        elif x > 1:
            Bo = Bo + L*numpy.dot(numpy.transpose(A),A)
            Wo = Wo + numpy.dot(numpy.transpose(C),C)
    
    B = (1.0/(len(tgrp)-1))*Bo
    W = (1.0/(X.shape[0] - len(tgrp)))*Wo
    
    return B,W
    
class MdAnova:
  def __init__(self):
    self.dimension = -1
    self.data = MdDatamatrix()
  def AddDataset(self,dataset):
    self.data.AddDataset( dataset )
    self.groupinfo_idx = 0
    self.group_list = {}
    i = 0
    for object in self.data.dataset.objects:
      group_key = object.group_list[self.groupinfo_idx]
      print group_key
      if( not self.group_list.has_key(group_key) ):
        self.group_list[group_key] = []
      self.group_list[group_key].append( i )
      i+=1
    return
  def Analyze(self):
    variable_index = 0
    data = {}
    data['sum'] = 0
    data['count'] = 0
    self.group_data = {}
    anova_data = []
    for key in self.group_list.keys():
      anova_data.append( [] )
      self.group_data[key] = {}
      within_group_sum = 0
      for observation_index in self.group_list[key]:
        val = self.data.matrix[variable_index,observation_index]
        anova_data[len(anova_data)-1].append( val )
        within_group_sum += val
        data['sum'] += val
      self.group_data[key]['sum'] = within_group_sum
      self.group_data[key]['count'] = len( self.group_list[key] )
      self.group_data[key]['average'] = ( self.group_data[key]['sum'] * 1.0 ) / ( 1.0 * self.group_data[key]['count'] )
      data['count'] += self.group_data[key]['count']
      print self.group_data[key]
    data['average'] = ( data['sum'] * 1.0 ) / ( data['count'] * 1.0 )
    bgss = 0
    wgss = 0
    for key in self.group_list.keys():
      bgss += self.group_data[key]['count'] * ( self.group_data[key]['average'] - data['average'] ) ** 2
      print key, bgss
      for observation_index in self.group_list[key]:
        val = self.data.matrix[variable_index,observation_index]
        wgss += ( val - self.group_data[key]['average'] ) ** 2
        print wgss
    df1 = ( len( self.group_list.keys() ) - 1 )
    bgms = bgss / df1
    print "bgms", bgms
    df2= ( data['count'] - len( self.group_list.keys() ) )
    wgms = wgss / df2
    print "wgms", wgms
    f = bgms/wgms
    print "F:", f
    prob = stats.fprob(df1,df2,f)
    print "p val", prob
    f, prob = stats.lF_oneway( anova_data[0], anova_data[1] )
    print f, prob
    return

class MdCanonicalVariate:
  def __init__(self):
    self.dimension = -1
    self.data = MdDatamatrix()
    
  def AddDataset(self,dataset):
    self.data.AddDataset( dataset )

  def Analyze(self):
    '''analyze'''
    #print "analyze"
    self.raw_eigen_values = []
    self.eigen_value_percentages = []
        
    #for d in self.datamatrix :
      #print d

    sums = []
    avrs = []
    ''' calculate the empirical mean '''
    for i in range ( self.data.nVariable ): 
      sums.append(  0 )
      for j in range ( self.data.nObservation ):
        sums[i] += self.data.matrix[i,j] 
    
    for sum in sums:
      avrs.append( float( sum ) / float( self.data.nObservation ) )
    
    #print "sum:", sums
    #print "avgs:",avrs
    #return
    
    for i in range ( self.nVariable ): 
      for j in range ( self.data.nObservation ):
        self.data.matrix[i,j] -= avrs[i] 

    #print self.datamatrix
    
    ''' covariance matrix '''
    self.covariance_matrix = numpy.dot( self.data.matrix, numpy.transpose(self.data.matrix)) / self.data.nObservation

    #print self.covariance_matrix
  
    ''' zz '''
    v, s, w = numpy.linalg.svd( self.covariance_matrix )
    #print "v", v
    #print "w", w
    
    #print "s[",
    self.raw_eigen_values = s
    sum = 0
    for ss in s:
      sum += ss
    for ss in s:
      self.eigen_value_percentages.append( ss/sum )
    cumul = 0
    eigen_values = []
    i = 0
    nSignificantEigenValue = -1
    nEigenValues = -1
    for ss in s:
      cumul += ss
      eigen_values.append( ss )
      #print sum, cumul, ss
      if cumul / sum > 0.95 and nSignificantEigenValue == -1:
        nSignificantEigenValue = i + 1
      if (ss /sum ) < 0.00001 and nEigenValues == -1:
        nEigenValues = i + 1
      i += 1
    
    #print nEigenValues, "eigen values obtained,", nSignificantEigenValue, "significant."
    #print eigen_values
    
    #for i in range( len(s) ):
      #print math.floor( ( s[i] / sum ) * 10000 + 0.5 ) / 100

    #print "s", int( s * 100 )/100
    #print "w", w
    #print v
    for i in range(nSignificantEigenValue):
      k = v[...,i]
      #print i, k, numpy.transpose(k) 
      det = numpy.dot( k, numpy.transpose( k ))
      #print det
    self.rotated_matrix = numpy.dot( w, self.data.matrix )
    #print w
    #print self.datamatrix[...,2]
    #print self.rotated_matrix[...,2]
    #print self.rotated_matrix
    self.loading = w
    
    self.new_dataset = self.data.dataset.copy()
    self.new_dataset.objects = []
    for i in range( self.nObservation ):
      object = MdObject()
      object.objname = self.data.dataset.objects[i].objname
      object.coords = self.rotated_matrix[...,i]
      object.group_list[:] = self.data.dataset.objects[i].group_list[:]
      self.new_dataset.objects.append( object )
      #if( i == 2 ) : print object.coords

class MdPrincipalComponent:
  def __init__(self):
    self.dimension = -1
    self.data = MdDatamatrix()
    #self.datamatrix = []
    return
  def AddDataset(self,dataset):
    self.data.AddDataset( dataset )

  def Analyze(self):
    '''analyze'''
    #print "analyze"
    self.raw_eigen_values = []
    self.eigen_value_percentages = []
        
    #for d in self.datamatrix :
      #print d

    sums = []
    avrs = []
    ''' calculate the empirical mean '''
    for i in range ( self.data.nVariable ): 
      sums.append(  0 )
      for j in range ( self.data.nObservation ):
        sums[i] += self.data.matrix[i,j] 
    
    for sum in sums:
      avrs.append( float( sum ) / float( self.data.nObservation ) )
    
    #print "sum:", sums
    #print "avgs:",avrs
    #return
    
    for i in range ( self.data.nVariable ): 
      for j in range ( self.data.nObservation ):
        self.data.matrix[i,j] -= avrs[i] 

    #print self.datamatrix
    
    ''' covariance matrix '''
    self.covariance_matrix = numpy.dot( self.data.matrix, numpy.transpose(self.data.matrix)) / self.data.nObservation

    #print self.covariance_matrix
  
    ''' zz '''
    v, s, w = numpy.linalg.svd( self.covariance_matrix )
    #print "v", v
    #print "w", w
    
    #print "s[",
    self.raw_eigen_values = s
    sum = 0
    for ss in s:
      sum += ss
    for ss in s:
      self.eigen_value_percentages.append( ss/sum )
    cumul = 0
    eigen_values = []
    i = 0
    nSignificantEigenValue = -1
    nEigenValues = -1
    for ss in s:
      cumul += ss
      eigen_values.append( ss )
      #print sum, cumul, ss
      if cumul / sum > 0.95 and nSignificantEigenValue == -1:
        nSignificantEigenValue = i + 1
      if (ss /sum ) < 0.00001 and nEigenValues == -1:
        nEigenValues = i + 1
      i += 1
    
    #print nEigenValues, "eigen values obtained,", nSignificantEigenValue, "significant."
    #print eigen_values
    
    #for i in range( len(s) ):
      #print math.floor( ( s[i] / sum ) * 10000 + 0.5 ) / 100

    #print "s", int( s * 100 )/100
    #print "w", w
    #print v
    for i in range(nSignificantEigenValue):
      k = v[...,i]
      #print i, k, numpy.transpose(k) 
      det = numpy.dot( k, numpy.transpose( k ))
      #print det
    self.rotated_matrix = numpy.dot( w, self.data.matrix )
    #print w
    #print self.datamatrix[...,2]
    #print self.rotated_matrix[...,2]
    #print self.rotated_matrix
    self.loading = w
    
    self.new_dataset = self.data.dataset.copy()
    self.new_dataset.objects = []
    for i in range( self.data.nObservation ):
      object = MdObject()
      object.objname = self.data.dataset.objects[i].objname
      object.coords = self.rotated_matrix[...,i]
      object.group_list[:] = self.data.dataset.objects[i].group_list[:]
      self.new_dataset.objects.append( object )
      #if( i == 2 ) : print object.coords
  